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Abstract 

When one deals with data drawn from continuous variables, a histogram is often inadequate to display their 
probability density. It deals inefhciently with statistical noise, and binsizes are free parameters. In contrast to that, 
the empirical cumulative distribution function (obtained after sorting the data) is parameter free. But it is a step 
fimction, so that its differentiation does not give a smooth probability density. Based on Fourier series expansion 
and Kolmogorov tests, we introduce a simple method, which overcomes this problem. Error bars on the estimated 
probability density are calculated using a jackknife method. We give several examples and provide computer code 
reproducing them. You may want to look at the corresponding figures 4 to 9 first. 
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1. Introduction 

One is often confronted with displaying an em- 
pirical probability density (PD) f{x) of a contin- 
uous variable x. Most commonly this is done us- 
ing histograms. While that is entirely appropriate 
and parameter free when x is confined to discrete 
values, in case of a continuous variable x, one is 
faced with choosing a binsize, or even several bin- 
sizes. This is a frustrated problem: On one side one 
would like the binsize infinitesimally small, so that 
the resolution of the curve becomes perfect. On 
the other side the binsize has to be sufficiently big 
so that statistical fluctuations do not destroy the 
smoothness of the underlying curve. Here we do 
not attempt to fine-tune the binsize parameter (s), 
but propose to by-pass the entire problem by us- 
ing a method based on the cumulative distribution 



function (CDF) 

X 

F{x) = j dx'fix') . (1) 

— oo 

F{x) is a monotonically increasing function as f{x) 
cannot be negative. In a range with f{x) > 0, F(x) 
is strictly monotone. 

Given a time series of n real numbers (data), a 
parameter free estimate of the CDF, called empir- 
ical CDF (ECDF), is well-known [1,2]: The step 
function F{x) which, after sorting the data, in- 
creases by 1/n at each data point. Unfortunately, 
that does not help directly in getting a good esti- 
mate of the probability density. The derivative of 
a step function is a sum of Dirac delta functions, 
whereas the probability density is often known to 
be a smooth function, as will be assumed in the 
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forthcoming. So it appears that one needs some 
kind of interpolation of the CDF before taking the 
derivative. This is no fun, as one has to decide 
whether the interpolation of 2, 3, 4, or k points will 
work best. In contrast, plotting a histogram of the 
data is simple and robust. However, it is tedious to 
guess a smooth function from a histogram. Typi- 
cally, cither the statistical errors on the bins are 
small, but the resolution of the function is bad, or 
the resolution is good, but the statistical errors are 
large. 

Let Fo{x) be the straight line, which is zero for 
X < Xtti and one for x > .t^^ , where is the 
smallest and .T^r^ the largest data point. After sub- 
tracting Fo{x) from F{x), a fmiction is left over, 
which is zero for x < and x > a;^„, and in- 
between well-suited for a Fourier series expansion. 
This leads to the desired smooth approximation 
as long as the expansion is sufficiently short, but 
will imitate every wiggle of the data, when carried 
too far. Therefore, one needs a cut-off criterion. 
We provide this by using the Kolmogorov [3,1,2] 
test, which tells us whether the difference between 
the ECDF and an analytical approximation of the 
CDF is explained by chance. The result is a well- 
defined, smooth empirical estimate of the PD. Our 
approach is so simple and straightforward, that one 
can hardly imagine that it is original, but we have 
not seen it in use before. Whether there exists pre- 
vious literature on it or not, it is certainly desirable 
to bring it to the attention of a general science com- 
munity. With this paper we also distribute startup 
software, which should help to bridge initial barri- 
ers against applying the approach. 

In the next section we review the CDF, its em- 
pirical estimate, and other preliminaries. Section 3 
explains our construction of PDs and gives nu- 
merical examples. We consider independent events 
sampled from (1) a Gaussian distributions, (2) a 
Cauchy distribution, and data from (3) a Markov 
chain Monte Carlo simulation, which creates an 
autocorrelated time series. Summary and conclu- 
sions follow in section 4. The appendix gives a list- 
ing of the main routine and explains Web access to 
Fortran 77 code, which reproduces the examples of 
this paper. 
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Fig. 1. ECDF from 100 Gaussian random numbers versus 
exact Gaussian CDF. 

2. Cumulative Distribution Functions 

Assume we generate n random numbers xi, ■ ■ ■, 
Xn according to a distribution function F{x). We 
may re-arrange the Xi in increasing order. Denoting 
the smallest value by a;^^ , the next smallest by ; 
etc., we arrive at 

X-jTi ^ ^7r2 — • • • — X-jr^ (2) 

where tti, . . . , 7r„ is a permutation of 1, . . . , n. As 
long as the data arc not yet created, the .x^^ are 
random variables, afterwards they are data. An es- 
timator for the distribution function F{x) is the 
ECDF 

^ 

F{x) = - for x„. <x < x^, , 1 , (3) 
n 

with i = 0, 1, . . . , n— 1, n and the definitions Xjro = 
—00, X7r„+i = +0O. Confidence limits can be ob- 
tained from the bimodal distribution: With Xq de- 
fined by F{xq) = q, < q < 1, the probability to 
find k data points with values < ■ . ■ < x.^^ < 
Xq and n — k data points with Xq < x,rfc+i < • • • < 
x,r„ is given by 

bn{k,q)=(J)qHl-qr-'' ■ (4) 

Figure 1 shows an ECDF from 100 Gaussian dis- 
tributed random numbers together with the exact 
CDF. Using Marsaglia [4] random numbers, the 
data were created for the probability density 
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Fig. 2. Peaked ECDF from 100 Gaussian random numbers 
versus exact Gaussian peaked CDF. 



cxp 



(5) 



and sorted with the heapsort algorithm. The CDF 
is in this case determined by the error function: 



G{x) 



dx'q(x') = + ^ erf , ^ 
^ 2 2 V\/2 



- 



(6) 



The computer code used was that of Ref . [2] . 

In contrast to histograms as estimators of a PD 
from data, the ECDF has the advantage that no 
free parameters are involved in its definition, but 
the probability density of events is encoded in its 
slope. This makes it often impossible to read off 
from graphs like Fig. 1 high probability regions, in 
particular, the point of maximum likelihood. An 
ECDF is well-suited for determining confidence in- 
tervals, however, and that can be further improved 
by switching from the CDF to the peaked CDF [2] 
defined by 



Fp{x) 



F{x) for F{x) < \ ; 
1 - F{x) for F(x) > 



(7) 



By construction the maximum of the peaked CDF 
Fp{x) is at the median xi and Fp{xi/2) = 1/2. 
Therefore, Fp{x) has two advantages over the CDF 
F{x): The median is clearly exhibited and the ac- 
curacy of the ordinate is doubled. It looks a bit like 
a PD, but it is in essence still the integrated PD. 

The peaked ECDF Fp (./;) is defined in the same 
way, just replacing F{x) in Eq. (7) by F{x). Fig- 



Fig. 3. Peaked ECDF from 10 000 Gaussian random num- 
bers versus exa<;t Gaussian peaked CDF. The arrows indi- 
cate 70% and 95% confidence intervals. 

ure 2 displays the peaked ECDF and peaked CDF 
for the same data as used in Fig. 1. While in these 
two figures deviations of the estimates from the 
exact function due to statistical deviations are 
clearly visible, the picture changes dramatically 
when one increases the number of data by a factor 
of 100. Figure 3 shows the peaked ECDF from 
10 000 Gaussian random numbers together with 
the exact peaked CDF. A difference is no longer 
visible to the naked eye. Note that for large n use 
of a fast sorting algorithm is mandatory. For a 
heapsort the CPU time scales with n log2(n). For 
details see, e.g., Ref. [2]. 

Estimation of confidence intervals is also illus- 
trated in Fig. 3: Just pick the desired likelihood on 
the ordinate and follow the arrows to their termi- 
nation points. Call such a point Xy, then the prob- 
ability for X < Xy is given by the value on the or- 
dinate for arrows emerging from the left, and the 
probability x > Xy for arrows emerging from the 
right. Using one arrow from the left, and one ar- 
row from the right, the probability content for the 
interval between these x values is obtained by sub- 
tracting the appropriate numbers from one (70% 
for the x values from the two inner and 95% for the 
X values from the two outer arrows of Fig. 3). 

Do the empirical and exact CDFs of our fig- 
ures agree? The Kolmogorov test [3] answers this 
question. It returns the probability Q, that the 
difference between the analytical CDF and an 
ECDF from statistically independent data is due 
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to chance. If the analytical CDF is exactly known 
and the data are indeed sampled from this dis- 
tribution, Q is a uniformly distributed random 
variable in the range < Q <1. Turned around, if 
one is not sure about the exact CDF, or the data, 
or both, and Q is small (say, Q < 10^®) one con- 
cludes that the difference between the proposed 
CDF and the data is not due to chance, but that 
one of the assumptions made is false. 

Kolmogorov's ingenious test relics on no more 
than the maximum difference A between the 
ECDF and the CDF. For the one-sided tests the 
relationship to Q is analytically known [5,2] for 
any value of n. In practicx^ (me prefers the two- 
sided test for which A is defined by 



A = max|F(2;) - F{x)\ 



(8) 



An exact analytical conversion into Q is then 
not known, but an asymptotic expansion due to 
Stephens [6] gives satisfactory results for n > 4, 
which is for practically all applications sufficient. 
The test yields Q = 0.19 for the data of Fig. 1 
and 2, and Q = 0.78 for the data of Fig. 3. Both 
values are consistent with the assumption that the 
difference between the data and the exact CDF is 
due to chance. 



3. Probability Densities 

We would like to construct an empirical prob- 
ability density (EPD) from an ECDF given by 
Eq. (3). The first idea, that comes to mind, is to 
differentiate smooth interpolations of F{x). This 
turns out to be tedious as long as one is unable to 
find a simple, generic rule, which determines the 
optimal interpolation for the purpose at hand, and 
we do not pursue this line. 

The method we propose consists of two well- 
defined steps. 

(i) Define as an initial approximation to F{x) 
a differcntiable, monotonically increasing 
function Fo{x), with details as specified be- 
low. 

(ii) Fourier expand the remainder until the Kol- 
mogorov test yields Q > Qcut = 1/2 (there 
may be some flexibility in lowering Qcut ) ■ 



For Fo{x) we require 

Fo{x) = for X < a, 
Fo{x) = 1 for a; > 6 , 



(9) 
(10) 



where the interval [a, b] has to lie within the range 
of the data: 



< a < b< x„ 



(11) 



For PDs with support on a compact interval, or 
with fast fall-off like for a Gaussian distribution, 
the natural choice is a = and b = Xt^^ . In case 
of slow fall-off, like for a Cauchy distributions, or 
other distributions with outliers, one has to restrict 
the analysis to [a, b] regions, which are sufficiently 
populated by data. This can be interpreted as con- 
sidering instead of /(x) the PD 



fab = I 



c/(.t) for a < a; < 6; 
otherwise . 



(12) 



Here the constant c is defined by the normaliza- 
tion / dx fab{x) = 1 and empirically obtained by 
the left-out probability content of F{x), i.e., c = 
n/riab, where Uab is the number of data in [a, b] and 
n the total number of data. We denote the CDF 
of fab{x) by Fab{x) aiid its ECDF by Fab{x). Af- 
ter calculating from Fab{x) an estimate fab{x) of 
fab{x), the estimate of f{x) is for x G [a,b] given 

hyJix)^c-'7M- 

For a: e [a, 6] we restrict our choice of Fo{x) in 
this paper to the straight line, 



X — Qi 

Fo{x) = for a<x<b. 

b — a 



(13) 



More elaborate definitions will likely give improve- 
ments in a number of situations (we comment on 
that in the conclusions), but would at the present 
state just discourage applications. Our point here 
is to show that good results are already obtained 
with the definition (13). 

Once Fo{x) is defined, the remainder of the 
ECDF is given by 

R{x)=Fab{x)-Foix), (14) 
which we expand into the Fourier series 

Rix) = ±d{i) sin (^^^^^) . (15) 

i=l ^ ^ 
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The cosine terms are not present due to the bound- 
ary conditions R{a) = R{b) = 0. The Fourier coef- 
ficients follow from 



d{i) 



dx R{x) 



sm 



ii7 (x — a) 



(16) 



Because in our case R{x) is the difference of a step 
function and a linear function, the integrals over 
the flat regions of the step function are easily cal- 
culated, and the integrals (16) are solved by adding 
them up. 

The Fourier expansion (14) is useless for a too 
large value of m, because it will then reproduce 
all statistical fluctuations of the data. To get 
around this problem, we perform the two-sided 
Kolmogorov test first between Fab{x) and Fq{x) 
(m = 0), and then each time m is incremented by 
m — > m + 1. Once Q > Qcut = 1/2 is reached for 
the Kolmogorov Q, we know that the difference 
between the data and our analytical approxima- 
tion is explained by statistical fluctuations. The 
other way round, the only information left in the 
data is statistical noise. Therefore, the expansion 
is terminated at that point. The thus obtained 
smooth estimate of the CDF, 



-F'estimate(a;) = Fo{x) + R{x) 



(17) 



yields fabi^') by differentiation, and our final esti- 
mate of the desired PD is: f(x) = Uabf abi^) l"^- 

One likes to attach error bars to the estimate 
of the PD. We do this by dividing the (unsorted) 
original data into jackknife [7,8] bins and repeat 
the analysis for each bin. Comparing the function 
values thus obtained at selected points, error bars 
follow in the usual jackknife way (see [2] for tech- 
nical details). 

3.1. Gaussian distribution 

Figure 4 shows a histogram of 51 bins for 2 000 
random numbers generated according to the Gaus- 
sian distribution (the error bars follow follow from 
the variance p ( 1 — p) of the bimodal distribution (4) 
with p = h{i)/n). 

Figure 5 shows the estimate g{x) obtained from 
the same data with the method described in this 



0.5 Histogram 
0.4 - 

^ 0.3 - : : 

0.2 - , H 

0.1 - 

L ' ' ' ' ' ' — 

-2 -1.5 -1 -0.5 0.5 1 1.5 2 

X 

Fig. 4. Histogram from 2 000 Gaussian random numbers. 
0.5 




Fig. 5. Estimate gix) of the Gaussian PD from the same 
data as used in Fig. 4. 

section. We used a = and b = a;,r„- For the 
estimate from all 2 000 data, Q = 0.97 was reached 
with m = 4 (Q = 0.056 with m = 3). Twenty 
jackknife bins were used to calculate the error bars. 
As all results of this section, the analysis is fully 
reproducible with the programs provided on the 
Web as described in the appendix. 

3.2. Cauchy distribution 

We consider the Cauchy distribution defined by 
the PD 



fcix) = 



7r(l-Fa;2) 
which leads to the CDF 



(18) 
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Fig. 6. Histogram from 20 000 Cauchy distributed random 
numbers. 




Fig. 7. Estimate /^(a;) of the Cauchy PD from the same 
data as used in Fig. 6. 

X 

Fc{x) = [ fc{x') dx' = l + - tan-1 {x) . (19) 

J 2 TT 



Due to the slow --^ x~'^ fall-off of the Cauchy PD, 
neither its mean nor its variance are defined. 

Therefore, it comes to no surprise that we en- 
counter considerably more difficulties in applying 
our method for Cauchy data than for Gaussian 
data or for the data of our next example. To over- 
come instabilities, the analysis has to be restricted 
to the central region and far more data than before 
are needed for stable estimates. We ended up with 
using 20 000 Cauchy distributed random number 
and restricting the analysis to 



a = X 



7T3001 ^^d b »^7ri7000 ' 



which came out to be a = —1.984 and b = 1.916. 
Figure 6 depicts the 51-bins histogram of the data 
in this region. The estimate of the PD /^{x) with 
our method is then shown in Fig. 7. As before, 
twenty jackknife bins were used for the error esti- 
mates and for the estimate from all data Q = 0.77 
was obtained for m = 2 (Q < 10~^^ for m = 1). 

3.3. Double peak from Markov Chain Monte Carlo 

Markov chain Monte Carlo (MCMC) simula- 
tions are widely employed in physics and other 
disciplines. Data in a MCMC time series are auto- 
correlated, which makes straightforward applica- 
tion of the Kolmogorov test questionable. Here we 
illustrate for an example from lattice gauge theory 
(LGT) a way to deal with this problem. 

In 4D U(l) LGT, a double peak in the action 
has been observed on symmetric lattices [9] . This is 
characteristic for a first order phase transition, al- 
though the situation in 4D U(l) LGT is somewhat 
questionable due to the weakness of the transition 
and other circumstances. 

For this paper we have used the biased 
Metropolis-Heatbath algorithm of Ref. [10] to 
generate U(l) data at /3 = 1.007 on an 8"^ lattice. 
These are parameters appropriate for producing 
a double peak. In LGT and statistical physics a 
standard unit for the time series is one "sweep", 
which corresponds for sequential updating (as 
used in our simulations) to updating each dynam- 
ical variable once. We have generated a time series 
of 2 560 000 sweeps, which takes about twenty 
hours on a 2GHz PC. 

We calculated the integrated autocorrelation 
time Tint of our time series and found Tint ~ 1 280 
sweeps. As the effective number of statistically 
independent data in an autocorrelated time se- 
ries is given by [2] n/Tint, our measurement of 
Tint implies that we have generated approximately 
2 560 000/1200 = 2 000 independent measure- 
ments. In the spirit of data reduction, we then 
selected 2 000 action values, separated by steps of 
1 200 sweeps, from the time series. 

For these 2 000 measurements, the histogram 
with 51 entries is plotted in Fig. 8, where s is the 
action density. Figure 9 shows the estimate f{s) 
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Fig. 8. Histogram from 2 000 effectively independent ciction 
measurements in U(l) lattice gauge theory. 




Fig. 9. Estimate f{s) of the PD from the same data as 
used in Fig. 8. 

obtained from the same data with our method, 
using a = x^^^ and b = a:,r„ • For the estimate from 
all 2 000 data Q = 0.93 was reached with m = 5 
{Q = 0.026 with m = 4), and twenty jackknife 
bins were used to calculate the error bars. The 
improvement from Fig. 8 to Fig. 9 is as good as 
the corresponding improvement for the Gaussian 
distribution. 

This result sheds new light on an effect, which is 
well-known to workers in MCMC simulations, but 
has to our knowledge not been satisfactorily ex- 
plained, the amazing smoothing, which one obtains 
when one includes all autocorrelated data in the 
histogram. In our case this is 640 000 data points, 
as we have taken measurements every 4 sweeps. 



The histogram error bars then fall almost on top 
of the jackknife error bars of Fig. 9. How can that 
be when there is little or no additional informa- 
tion in the autocorrelated data? Part of the answer 
appears to be that most of this information is al- 
ready in the statistically independent data, but is 
usually not exploited. 



4. Summary and Conclusions 

To someone used to plotting histogram, our 
method may appear difficult, but indeed it is not. 
True, first one has to sort the data to calculate 
their ECDF, then from it the empirical PD, using 
Fourier scries expansion and Kolmogorov tests, 
and finally jackknife error bars need to be calcu- 
lated, but all these steps are standard. Besides 
the subroutine shown in the appendix, another 
short subroutine, and the main program, we had 
no programming to do. All other routines were 
already available in the Web based Fortran 77 
code of Ref. [2] . Once the archive described in the 
appendix is downloaded from the Web, and the 
program is up and running, the extra effort com- 
pared to plotting histograms is almost negligible. 
And the differences between Fig. 4 and Fig. 5, 
Fig. 6 and Fig. 7, and finally Fig. 8 and Fig. 9 
speak for themselves. 

It is astonishing that these results could be ob- 
tained with a simple straight line (13) as initial 
approximation for the CDF. There is certainly 
space for improvement at the price of giving up 
the presently achieved simplicity. For our exam- 
ples one and two, the exact CDFs are known and 
pointed out by Eqs. (6) and (19). Obviously, one 
does not test our method by constructing ^0(2^) 
from them in these examples, but they could be 
good choices, when one has reasons to expect them 
to be a close, albeit not exact, approximation of 
the investigated distribution. For instance, for the 
double peak of our third example, a good initial 
guess could be a sum of two Gaussians, and one 
could start off by fitting their CDF to the data. 
Afterwards care has to be taken that Eqs. (9) 
and (10) remain valid. That this can always be 
achieved follows from the interpretation (12) of 
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that requirement. We abstained from investigating 
a double Gaussian Fo{x), because it starts to get 
tedious and there is nothing really to improve on 
our result. But, one can easily imagine that there 
are more complicated situations, accompanied by 
limited statistics, where an improvement of Fq(x) 
becomes essential. If the problem for which this 
happens is important too, it will become worth- 
while to explore more Fq{x) functions. 

With our Qcut = 1/2 rule, we are slightly over- 
expanding the Fourier transformation (15). In the 
average Q should be 1/2, but all our final values are 
> 1/2. That gives some flexibility to lower Qcut, 
which should be used with discretion in situations 
were the m of the Fourier expansion (15) appears 
to be too large. A warning right away: The only 
situation in which we did not see a rapid approach 
towards Q > Qcut turned out to be one in which 
we had not noticed that the underlying distribu- 
tion was discrete, while the Kolmogorov test did 
notice that. 

We did not develop a statistically rigorous ap- 
proach. We address physicists and others, who do 
not hesitate to use whatever works, not those, who 
want to forbid numerical recipes. We rely on the 
assumption that the Fourier series (15) would be 
rapidly convergent, when the ECDF in Eq. (14) 
would be replaced by the corresponding (unknown) 
exact CDF. That is in spirit similar to picking a 
primer in Bayesian statistics, when a rigorous one 
is not known, a procedure, which can modify the 
probability content of confidence intervals. 
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Appendix A. Computer Code 

An archive with Fortran 77 example runs can 
presently be downloaded from the website of BB. 
Google Bernd Berg, or go directly to 



http : //people/scs . f su . edu/~ berg/ . 



Take from there the research link and download 
the gzipped archive 



STMC.CDFtoPD.tgz . 

Unfold the archive. (If instructions for that are 
needed, they can be found on the website of Ref . [2] , 
which is also linked on the main website of BB.) 
Go then to the folder 

STMC_CDFtoPD/Work/CDFtoPDexamples . 

All three examples of this paper can be run as 
special cases of the program cdf _to_pd. f in this 
folder, and more instructions are given in its 
readme.txt file. To show that at its heart our 
method is indeed quite simple, we list in the fol- 
lowing our main subroutine. 

SUBROUTINE CDF_PD (lUD ,NDAT , SDAT , Fxct ) 
C Bernd Berg, Robert Harris Dec 16 2007. 

C Transforms an empirical cumulative distribution function (CDF) into 
C a corresponding probability density using and initial function plus 
C Fourier series expansion for the CDF, 
C On INPUT: 

C lUO Write unit , unchanged on exit , 
C NDAT Number of input data, unchanged on exit. 
C SDAT Sorted input data, unchanged on exit. 
C INTERNAL: 

C Qcut Cut-off, Fourier expansion is terminated for Q>Qcut , 

C cumulative distribution function, 

C NMAX Maximum number of terms in the Fourier series, 

C On OUTPUT: 

C Fxct Exact CDF (means here analytical approximation of the CDF) , 
include ' , ,/, , /Libs/Fortran/implicit , sta' 
include ' , ,/, , /Libs/Fortran/constants ,par ' 
PARAMETERCQcut-HALF,NMAX=100) ! Change also in functions, 
DIMENSION SDAT(NDAT) ,Fxct(NDAT) 

COMMON /CDFProb/ XMIN,XRANGE,DN(NMAX) ,M ! Expansion parameters, 

C 

XHIN-SDATCl) ! Initializations, 
XRANGE=SDAT (NDAT) -SDAT ( 1 ) 
DO J=1,NDAT 

Fxct(J)-(SDAT(J)-SDAT(l))/XRANGE 
END DO 
DO K-1,NHAX 

DN(K)=ZERO 
ENDDO 

DO M=1,NMAX ! Integration for the Fourier series coefficients: 
DO I-1,NDAT-1 

X1-(SDAT(1)-SDAT(1))/XRANGE 
X2-(SDAT(1+1)-SDAT(1)) /XRANGE 

DN(M)=DN(M)+(DNE»I/NDAT-X1)*C0SCM»P1»X1)/(M*P1) 
DN(M)-DN(M) - (DNE*I/NDAT-X2) *COS (M*P1»X2) / (M*P1) 
DNCM)-DNCM)+S1NCM»P1*X1)/(M*M»P1*PI) 
DN(M)-DN(M)-S1N(H*P1*X2)/(M*M*P1*PI) 
ENDDO 

DN(M)-DN(M)*TWO 

DO K=1,NDAT ! CDF in Fourier series approximations; 
XRANGE1=SDAT (K) -SDAT (1 ) 

Fxct(K)-Fxct (K)+DNCM)»S1N(M*PI/XRANGE»XRANGE1) 
ENDDO 

CALL KDLM2_ASCNDAT, Fxct, DEL, 0) ! Kolmogorov test, 
IF(Q,GT,Qcut) GOTO 1 
ENDDO 

WRITE(IUO,'(/," CDF.PD failed M,q =" , I6,G12 , 3) ' ) M.Q 

STOP "CDF.PD: Expansion failed," 
1 WRITE(IUO,'(/," CDF.PD: Final M,q -" , 16 ,G12 , 3, /) " ) M.Q 

C 

RETURN 
END 

With exception of the jackknifc routine 



dat_to_dat j .f . 
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to be found in Libs/Fortran of the package, all 
other subroutines used are from the code of Ref . [2] 
and for convenience included here. For the Kol- 
mogorov test in the form of Stephens the routine 
kolm2_as . f of [2] has been modified to abort in 
case of no convergence. (Note also that the related 
routine kolm2_as2.f of [2], which is not used in 
this paper, does not work for large values of the in- 
put arguments Ni , A'2 due to bad coding of a mul- 
tiplication of Ni and N2, which can be easily cor- 
rected.) All routines arc copyrighted by their au- 
thors, who are listed in one of the first lines of each 
routine. Limited permission of their use is given 
under the conditions stated on the Fortran down- 
load page of Ref. [2] . 
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